#include<iostream>
#include<limits>
#include<string>
#include<stdlib.h>
#include<cmath>
using namespace std;

#define M 50
const double epsl=100*numeric_limits<double>::epsilon();

int perm(int k,int n)
{
   int i,m=1;
   for(i=n;i>=1;i--)
   {
      if(i>=k+1)
          m=m*i;
      if(i<=n-k)
          m=m/i;
   }
   return m;
}

class Function
{
   public:
           virtual double operator()(const double &x)const=0;
           virtual double diff(const int &order,const double &x)const
           {
               double d=0;
               int i;
               for(i=0;i<=order;i++)
               {
                   d=d+perm(i,order)*pow(-1,i)*(*this)(x+(order-2*i)*epsl);
               }
               return d/pow(2*epsl,order);
           }
};

class Subroutine
{
   private:
           Function &f;
           Function &g;
           double pro,sum;
           int i,j,n,k;
           double I[M];
   public:
           double c[M];
           Subroutine(Function &f,Function &g,int n):f(f),g(g),n(n){}
           void generate()
           {
               for(k=0;k<=n;k++)
               {
                   I[k]=g(double(k));
               }              
           }
           double Newtvalue(double x)
           {
               double newt[n+1][n+1];
               for(i=0;i<=n;i++)
               {
                   newt[i][0]=f(I[i]);
               }
               pro=1;
               sum=newt[0][0];
               c[0]=newt[0][0];
               for(i=1;i<=n;i++)
               { 
                   for(j=i;j<=n;j++)
                   {
                       newt[j][i]=(newt[j][i-1]-newt[j-1][i-1])/(I[j]-I[j-i]);
                   }
                   pro=pro*(x-I[i-1]);
                   sum=sum+pro*newt[i][i];
                   c[i]=newt[i][i];
               }
               return sum;
            }
            double Hervalue(double x)
            {
               double her[n+1][n+1]={0};
               int Order[n+1]={0};
               for(i=1;i<=n;i++)
               {
                   if(I[i]==I[i-1])
                       Order[i]=Order[i-1]+1; 
               }
               for(i=0;i<=n;i++)
               {
                   for(j=0;j<=Order[i];j++)
                   {
                       if(j==0)
                           her[i][j]=f(I[i]);
                       else
                           her[i][j]=f.diff(j,I[i]);
                   }
               }
               pro=1;
               sum=her[0][0];
               c[0]=her[0][0];
               for(i=1;i<=n;i++)
               { 
                   for(j=i;j<=n;j++)
                   {
                       if(her[j][i]==0)
                           her[j][i]=(her[j][i-1]-her[j-1][i-1])/(I[j]-I[j-i]);
                   }
                   pro=pro*(x-I[i-1]);
                   sum=sum+pro*her[i][i];
                   c[i]=her[i][i];
               }
               return sum;
            }
};
